import pandas as pd
import numpy as np
import scipy.interpolate as spi
import matplotlib.pyplot as plt
from smooth_spline import x_to_arc
from rotate_and_interp import get_leading, leading_interpolation
from bspline_process import translation
from utils2 import set0

def judge_Monotonicity(ordinates):
    '''
        判断单调性（递增）
    :param ordinates:
    :return:
    '''

    for i in range(ordinates.shape[0]):
        if ordinates[i, 0] > ordinates[i+1, 0]:
            print(ordinates[i:i+10])
            return False

    return True

def produce_x_v1(N=65, keep_decimals=False):
    '''
        产生余弦分布的x序列(上表面65个点；下表面不包括前缘顶点，为64个点)
    :return:
    '''
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    phi = np.pi / 2 / (N - 1)

    x = 1 - np.cos(np.arange(0, N, 1) * phi)
    reverse_x = x[::-1]

    res = np.concatenate((reverse_x, x[1::]), axis=0)

    return res


def produce_x_v2(N=65, keep_decimals=False):
    '''
        产生余弦分布的x序列(上表面65个点；下表面不包括前缘顶点，为64个点)
        为了适应改进插值方法，要将横坐标转换成弧长的形式
    :return:
    '''
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    phi = np.pi / 2 / (N - 1)

    x = 1 - np.cos(np.arange(0, N, 1) * phi)
    reverse_x = x[::-1]

    res = np.concatenate((reverse_x, x[1::]), axis=0)

    return res


def produce_x_v3(N=65, keep_decimals=False):
    '''
        产生余弦分布的x序列(上表面65个点；下表面不包括前缘顶点，为64个点)
        0.9999...置为1
    :return:
    '''
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    phi = np.pi / 2 / (N - 1)

    x = 1 - np.cos(np.arange(0, N, 1) * phi)
    x[-1] = 1
    reverse_x = x[::-1] * -1

    res = np.concatenate((reverse_x, x[1::]), axis=0)

    return res


def produce_x_v4(N=65, keep_decimals=False):
    '''
        产生余弦分布的x序列(上表面65个点；下表面不包括前缘顶点，为64个点)
        0.5附近较密集
        0.9999...置为1
    :return:
    '''
    # 保留3位小数
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    phi = np.pi / 2 / (N - 1)

    x = 1 - np.cos(np.arange(0, N, 1) * phi)
    x[-1] = 1

    # 为了得到一个在0.5附近较密集，其余部分稀疏的数据，将进行以下操作
    # 0-0.5部分由原分布乘上(-0.5)再加0.5
    # 0.5-1部分由原分布乘上(0.5)再加0.5

    x1 = x * (-0.5) + 0.5
    x2 = x * 0.5 + 0.5

    res = np.concatenate([x1, x2], axis=0)

    return res


def produce_x_v5(N=65, keep_decimals=False, isUp=True):
    '''
        产生余弦分布的x序列(上表面65个点；下表面不包括前缘顶点，为64个点)
        0.9999...置为1
    :return:
    '''
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    phi = np.pi / 2 / (N - 1)

    x = 1 - np.cos(np.arange(0, N, 1) * phi)
    x[-1] = 1
    reverse_x = x[::-1]

    if isUp:
        return reverse_x
    else:
        return x


def produce_x_v6(N=65, keep_decimals=False):
    '''
    :return:
    '''
    # 保留3位小数
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    tN = N  # 2 * N - 1
    x = np.power((1 - np.cos((np.arange(0, tN, 1)) * np.pi / (N - 1))), 2) / 4

    # 为了得到一个在0.5附近较密集，其余部分稀疏的数据，将进行以下操作
    # 0-0.5部分由原分布乘上(-0.5)再加0.5
    # 0.5-1部分由原分布乘上(0.5)再加0.5

    x1 = x * (-0.5) + 0.5
    x2 = x * 0.5 + 0.5

    res = np.concatenate([x1[::-1], x2], axis=0)

    return res


def produce_x_v7(N=65, keep_decimals=False):
    '''
    :return:
    '''
    # 保留3位小数
    if keep_decimals == True:
        np.set_printoptions(precision=3)

    tN = N  # 2 * N - 1
    x = np.power((1 - np.cos((np.arange(0, tN, 1)) * np.pi / (N - 1))), 2) / 4

    return x


def produce_x_v8(N=65):
    '''
        产生余弦分布的x序列(上表面151(N)个点；下表面不包括前缘顶点，为150(N-1)个点。另外，上下表面都必须包括后缘点)
    :return:
    '''

    tN = 2 * N - 1
    x = np.power((1 - np.cos((np.arange(0, tN, 1) - (N - 1)) * np.pi / (N - 1))), 2) / 4

    return np.concatenate([x[0:N], x[N:2 * N - 1]], axis=0)


def cubic(start, end, *args):
    '''三次样条插值'''

    count = 4 * (len(args) - 1)
    mat_ori = np.zeros((count, count))
    mat_ans = np.zeros(count)
    index = 0
    for i, j in zip(args[:-2], args[1:-1]):
        mat_ori[4 * index, 4 * index] = i[0] ** 3
        mat_ori[4 * index, 4 * index + 1] = i[0] ** 2
        mat_ori[4 * index, 4 * index + 2] = i[0]
        mat_ori[4 * index, 4 * index + 3] = 1
        mat_ans[4 * index] = i[1]
        mat_ori[4 * index + 1, 4 * index] = j[0] ** 3
        mat_ori[4 * index + 1, 4 * index + 1] = j[0] ** 2
        mat_ori[4 * index + 1, 4 * index + 2] = j[0]
        mat_ori[4 * index + 1, 4 * index + 3] = 1
        mat_ans[4 * index + 1] = j[1]
        mat_ori[4 * index + 2, 4 * index] = 3 * j[0] ** 2
        mat_ori[4 * index + 2, 4 * index + 1] = 2 * j[0]
        mat_ori[4 * index + 2, 4 * index + 2] = 1
        mat_ori[4 * index + 2, 4 * index + 4] = -3 * j[0] ** 2
        mat_ori[4 * index + 2, 4 * index + 5] = -2 * j[0]
        mat_ori[4 * index + 2, 4 * index + 6] = -1
        mat_ans[4 * index + 2] = 0
        mat_ori[4 * index + 3, 4 * index] = 6 * j[0]
        mat_ori[4 * index + 3, 4 * index + 1] = 2
        mat_ori[4 * index + 3, 4 * index + 4] = -6 * j[0]
        mat_ori[4 * index + 3, 4 * index + 5] = -2
        mat_ans[4 * index + 3] = 0
        index += 1
    mat_ori[4 * index, 4 * index] = args[-2][0] ** 3
    mat_ori[4 * index, 4 * index + 1] = args[-2][0] ** 2
    mat_ori[4 * index, 4 * index + 2] = args[-2][0]
    mat_ori[4 * index, 4 * index + 3] = 1
    mat_ans[4 * index] = args[-2][1]
    mat_ori[4 * index + 1, 4 * index] = args[-1][0] ** 3
    mat_ori[4 * index + 1, 4 * index + 1] = args[-1][0] ** 2
    mat_ori[4 * index + 1, 4 * index + 2] = args[-1][0]
    mat_ori[4 * index + 1, 4 * index + 3] = 1
    mat_ans[4 * index + 1] = args[-1][1]
    mat_ori[4 * index + 2, 0] = 3 * args[0][0] ** 2
    mat_ori[4 * index + 2, 1] = 2 * args[0][0]
    mat_ori[4 * index + 2, 2] = 1
    mat_ans[4 * index + 2] = start
    mat_ori[4 * index + 3, 4 * index] = 3 * args[-1][0] ** 2
    mat_ori[4 * index + 3, 4 * index + 1] = 2 * args[-1][0]
    mat_ori[4 * index + 3, 4 * index + 2] = 1
    mat_ans[4 * index + 3] = end
    mat_rg = np.linalg.solve(mat_ori, mat_ans)

    def rtn_func(x):
        def bin_search(left, right):
            if (x < args[left][0]): return left
            if (x > args[right - 1][0]): return right - 1
            if (right - left <= 1): return left
            mid = int((left + right) / 2)
            if (x < args[mid][0]):
                return bin_search(left, mid)
            else:
                return bin_search(mid, right)

        num = bin_search(0, len(args) - 1)
        return mat_rg[4 * num] * x ** 3 + mat_rg[4 * num + 1] * x ** 2 + mat_rg[4 * num + 2] * x + mat_rg[4 * num + 3]

    return rtn_func


def Cubic_Spline_Interpolation_v1(ordinates, x, N=65):
    '''
        做三次样条插值
    :param ordinates: 坐标集合
    :param x: 插值点
    :return: 插值后的坐标集合
    '''
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

    # 前缘顶点
    key = np.argmin(ordinates[:, 0])

    # 上表面坐标点
    ordinates_up = ordinates[0:key + 1]
    ordinates_up = ordinates_up[::-1]

    # 下表面坐标点
    ordinates_low = ordinates[key::]

    # 进行三次样条拟合
    # 上机翼表面
    # spi.interp1d
    # method1
    # print(ordinates_up)
    ipo3_u = spi.splrep(ordinates_up[:, 0], ordinates_up[:, 1], k=3)  # 样本点导入，生成参数
    iy3_u = spi.splev(x[0:N], ipo3_u)  # 根据观测点和样条参数，生成插值
    # method2
    # spi.CubicSpline()
    # f1 = spi.interp1d(ordinates_up[:, 0], ordinates_up[:, 1], kind='cubic')
    # iy3_u = f1(x[0:N + 1])

    # 下机翼表面
    # method1
    ipo3_l = spi.splrep(ordinates_low[:, 0], ordinates_low[:, 1], k=3)  # 样本点导入，生成参数
    iy3_l = spi.splev(x[N - 1::], ipo3_l)  # 根据观测点和样条参数，生成插值

    # method2
    # f2 = spi.interp1d(ordinates_low[:, 0], ordinates_low[:, 1], kind='cubic')
    # iy3_l = f2(x[N + 1::])

    # 新的坐标集合
    ordinates_interpoaltion = np.empty(shape=(2 * N - 1, 2))
    ordinates_interpoaltion[:, 0] = x
    ordinates_interpoaltion[0:N, 1] = iy3_u
    ordinates_interpoaltion[N::, 1] = iy3_l[1::]

    # visualize
    # plt.scatter(ordinates[0:key + 1, 0], ordinates[0:key + 1, 1], c='b', label="原机翼")
    # plt.scatter(ordinates[key + 1::, 0], ordinates[key + 1::, 1], c='b', label="原机翼")
    # plt.scatter(ordinates_interpoaltion[0:N + 1, 0], ordinates_interpoaltion[0:N + 1, 1], c='r', label="插值生成的机翼")
    # plt.scatter(ordinates_interpoaltion[N + 1::, 0], ordinates_interpoaltion[N + 1::, 1], c='r', label="插值生成的机翼")
    #
    # plt.ylim(-0.3, 0.3)
    # plt.ylabel('指数')
    # plt.title('机翼数据三次样条插值拟合结果')
    # plt.legend()
    # plt.show()

    return ordinates_interpoaltion


def Cubic_Spline_Interpolation_v2(ordinates, x, N=65):
    '''
        做三次样条插值
    :param ordinates: 坐标集合
    :param x: 插值点
    :return: 插值后的坐标集合
    '''
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

    # 前缘顶点
    key = np.argmin(ordinates[:, 0])

    # 上表面坐标点
    ordinates_up = ordinates[0:key + 1]
    ordinates_up = ordinates_up[::-1]

    # 下表面坐标点
    ordinates_low = ordinates[key::]

    # 横坐标转换为弧长
    arc_ordinates_up, arc_ordinates_low, arc_ordinates = x_to_arc(ordinates_up, ordinates_low, isShow=False)
    # print(arc_ordinates)
    arc_min = np.min(arc_ordinates[0:key + 1, 0])
    arc_max = np.max(arc_ordinates[key::, 0])

    # 对弧长进行三次样条拟合
    # x放缩到弧长的范围
    x[0:N] = x[0:N] * arc_min * -1
    x[N::] = x[N::] * arc_max
    # 由于机翼上下表面已经被“拉”成一条曲线，所以可以直接拿这条曲线作插值
    # spi.interp1d
    # 对x作插值
    ipo3_x = spi.splrep(arc_ordinates[:, 0], ordinates[::-1, 0], k=3)  # 样本点导入，生成参数
    iy3_x = spi.splev(x, ipo3_x)  # 根据观测点和样条参数，生成插值

    # 对y作插值
    ipo3_y = spi.splrep(arc_ordinates[:, 0], ordinates[::-1, 1], k=3)  # 样本点导入，生成参数
    iy3_y = spi.splev(x, ipo3_y)  # 根据观测点和样条参数，生成插值

    # 新的坐标集合
    ordinates_interpoaltion = np.empty(shape=(2 * N - 1, 2))
    ordinates_interpoaltion[:, 0] = iy3_x
    ordinates_interpoaltion[:, 1] = iy3_y
    # print(x.shape)
    # print(ordinates_interpoaltion)

    # visualize
    # plt.scatter(ordinates[:, 0], ordinates[:, 1], c='b', s=5, label="原机翼")
    # plt.scatter(ordinates_interpoaltion[:, 0], ordinates_interpoaltion[:, 1], c='r', s=5, label="插值生成的机翼")
    # plt.ylim(-0.3, 0.3)
    # plt.ylabel('指数')
    # plt.title('机翼数据三次样条插值拟合结果')
    # plt.legend()
    # plt.show()
    # plt.clf()

    # 弧长-x
    # plt.scatter(arc_ordinates[:, 0], ordinates[::-1, 0], c='b')
    # plt.scatter(x, iy3_x, c='r')
    # plt.show()
    # plt.clf()

    # 弧长-y
    # plt.scatter(arc_ordinates[:, 0], ordinates[::-1, 1], c='b')
    # plt.scatter(x, iy3_y, c='r')
    # plt.show()
    # plt.clf()

    return ordinates_interpoaltion


def Cubic_Spline_Interpolation_v3(ordinates, x, N=65):
    '''
        对(弧长,x)，(弧长,y)做三次样条插值
        存在的问题是，机翼前缘存在内勾的情况
    :param ordinates: 坐标集合
    :param x: 插值点
    :return: 插值后的坐标集合
    '''
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

    # 前缘顶点
    key = np.argmin(ordinates[:, 0])

    # 上表面坐标点
    ordinates_up = ordinates[0:key + 1]
    ordinates_up = ordinates_up[::-1]

    # 下表面坐标点
    ordinates_low = ordinates[key::]

    # 横坐标转换为弧长
    arc_ordinates_up, arc_ordinates_low, arc_ordinates = x_to_arc(ordinates_up, ordinates_low, isShow=False)
    # print(arc_ordinates)
    arc_min = np.min(arc_ordinates[0:key + 1, 0])
    arc_max = np.max(arc_ordinates[key::, 0])

    # 对弧长进行三次样条拟合
    # x放缩到弧长的范围
    x[0:N] = x[0:N] * arc_min * -1
    x[N::] = x[N::] * arc_max
    # 由于机翼上下表面已经被“拉”成一条曲线，所以可以直接拿这条曲线作插值
    # spi.interp1d
    # 对x作插值
    ipo3_x = spi.splrep(arc_ordinates[:, 0], ordinates[::-1, 0], k=3)  # 样本点导入，生成参数
    iy3_x = spi.splev(x, ipo3_x)  # 根据观测点和样条参数，生成插值

    # 对y作插值
    ipo3_y = spi.splrep(arc_ordinates[:, 0], ordinates[::-1, 1], k=3)  # 样本点导入，生成参数
    iy3_y = spi.splev(x, ipo3_y)  # 根据观测点和样条参数，生成插值

    # 新的坐标集合
    ordinates_interpoaltion = np.empty(shape=(2 * N - 1, 2))
    ordinates_interpoaltion[:, 0] = iy3_x
    ordinates_interpoaltion[:, 1] = iy3_y

    # visualize
    plt.scatter(ordinates[:, 0], ordinates[:, 1], c='b', s=5, label="原机翼")
    plt.scatter(ordinates_interpoaltion[:, 0], ordinates_interpoaltion[:, 1], c='r', s=5, label="插值生成的机翼")
    plt.ylim(-0.3, 0.3)
    plt.ylabel('指数')
    plt.title('机翼数据三次样条插值拟合结果')
    plt.legend()
    plt.show()
    plt.clf()

    # 弧长-x
    # plt.scatter(arc_ordinates[:, 0], ordinates[::-1, 0], c='b')
    # plt.scatter(x, iy3_x, c='r')
    # plt.show()
    # plt.clf()

    # 弧长-y
    # plt.scatter(arc_ordinates[:, 0], ordinates[::-1, 1], c='b')
    # plt.scatter(x, iy3_y, c='r')
    # plt.show()
    # plt.clf()

    return ordinates_interpoaltion


def B_Spline_Interpolation_v1(ordinates, fx, N=65):
    '''
        使用B样本进行插值（B样条插值点：0.5-0，0.5-1）
    :param ordinates: 原机翼坐标
    :param fx: 横坐标分布
    :param N: 上/下表面翼型点数
    :return:
    '''

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    np.set_printoptions(threshold=np.inf)  # 打印不省略

    # 定义N为上表面的数据点个数
    N = 65

    # get x, y
    x = ordinates[:, 0].flatten()
    y = ordinates[:, 1].flatten()

    # print("ordinates")
    # print(ordinates)
    # for it in ordinates:
    #     print('''[{},{}],'''.format(it[0], it[1]))

    # B样条插值
    tck, u = spi.splprep([x, y], s=0, per=0)  # 注：per==0不封闭，per==1封闭
    xi, yi = spi.splev(fx, tck)

    ordinates_bspline = np.empty(shape=(len(xi), 2))
    ordinates_bspline[:, 0] = xi
    ordinates_bspline[:, 1] = yi

    # 再次插值得到横坐标相同的分布
    # 插值方法：三次样条插值
    x_up = produce_x_v7(N=N)
    x_low = produce_x_v7(N=N)
    ordinates_CbSp = np.empty(shape=(N * 2 - 1, 2))

    # print("ordinates_bspline")
    # # print(ordinates_bspline)
    # for it in ordinates_bspline:
    #     print('''[{},{}],'''.format(it[0], it[1]))

    if ordinates_bspline[0, 1] > 0 or \
            len(np.where(ordinates_bspline == np.min(ordinates_bspline[:, 0]))[0]) == 2:
        # 保存为逆时针排列的翼型(最密集的地方在x轴上方)
        flag1 = 0
        for i in range(ordinates_bspline.shape[0] - 1):
            if ordinates_bspline[i, 0] > ordinates_bspline[i + 1, 0]:
                flag1 = i + 1
                break
        # flag2 = np.argmin(ordinates_bspline[:, 0]) # 只能找到第一个最小的
        flag2 = np.where(ordinates_bspline == np.min(ordinates_bspline[:, 0]))[0][-1]
        ordinates_bspline_tmp = np.concatenate([ordinates_bspline[0:flag1][::-1],
                                                ordinates_bspline[flag1 + 1:flag2 + 1],
                                                ordinates_bspline[flag2 + 1::]], axis=0)

    else:
        # 保存为逆时针排列的翼型(最密集的地方在x轴下方)
        # flag2 = np.argmin(ordinates_bspline[:, 0]) # 只能找到第一个最小的
        flag2 = np.where(ordinates_bspline == np.min(ordinates_bspline[:, 0]))[0][-1]

        flag1 = 0
        for i in range(flag2 + 1, ordinates_bspline.shape[0] - 1):
            if ordinates_bspline[i, 0] > ordinates_bspline[i + 1, 0]:
                flag1 = i + 1
                break

        ordinates_bspline_tmp = np.concatenate([ordinates_bspline[flag2:flag1][::-1],
                                                ordinates_bspline[0:flag2][::-1],
                                                ordinates_bspline[flag1 + 1::]], axis=0)

    # 将翼型规范化
    ordinates_bspline_tmp = translation(ordinates_bspline_tmp)

    # 翼型数据截断
    ordinates_bspline_tmp = set0(ordinates_bspline_tmp)
    # print("tmp")
    # print(ordinates_bspline_tmp)

    # 前缘顶点
    key = np.argmin(ordinates_bspline_tmp[:, 0])
    # print(ordinates_bspline_tmp)
    # 上表面坐标点
    ordinates_CbSp_up = ordinates_bspline_tmp[0:key + 1]
    ordinates_CbSp_up = ordinates_CbSp_up[::-1]

    # 下表面坐标点
    ordinates_CbSp_low = ordinates_bspline_tmp[key::]

    # pd.DataFrame(np.reshape(ordinates_CbSp_up, newshape=(-1, 2))).to_csv("test_bspline_up.csv", index=None, columns=None)
    # pd.DataFrame(np.reshape(ordinates_CbSp_low, newshape=(-1, 2))).to_csv("test_bspline_low.csv", index=None, columns=None)
    # pd.DataFrame(np.reshape(ordinates_bspline, newshape=(-1, 2))).to_csv("test_bspline.csv", index=None, columns=None)

    # print("low")
    # print(ordinates_CbSp_low)
    # print("up")
    # print(ordinates_CbSp_up)

    # 上表面
    ipo3_x = spi.splrep(ordinates_CbSp_up[:, 0], ordinates_CbSp_up[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[0:N, 1] = spi.splev(x_up, ipo3_x)  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[0:N, 0] = x_up

    # 下表面
    ipo3_x = spi.splrep(ordinates_CbSp_low[:, 0], ordinates_CbSp_low[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[N - 1::, 1] = spi.splev(x_low, ipo3_x)  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[N - 1::, 0] = x_low

    # 翼型数据截断
    ordinates_CbSp = set0(ordinates_CbSp)

    # show
    # plt.scatter(ordinates[:, 0], ordinates[:, 1], c='b', s=20, label="原机翼")
    # plt.scatter(ordinates_bspline[:, 0], ordinates_bspline[:, 1], c='r', s=5, label="B-样条拟合的机翼")
    # plt.scatter(ordinates_bspline_tmp[:, 0], ordinates_bspline_tmp[:, 1], c='yellow', s=5, label="B-样条拟合的机翼（规范化）")
    # plt.scatter(ordinates_CbSp[:, 0], ordinates_CbSp[:, 1], c='black', s=5, label="三次样条插值拟合的机翼")
    # plt.plot([0, 1], [0, 0], c='orange', linestyle="--")  # 横轴
    # plt.plot([0, 0], [-0.3, 0.3], c='orange', linestyle="--")  # 纵轴
    # plt.ylim(-0.3, 0.3)
    # plt.ylabel('指数')
    # plt.title('机翼数据B样条插值拟合结果')
    # plt.legend()
    # plt.show()

    return ordinates_CbSp


def B_Spline_Interpolation_v2(ordinates, fx, N=65):
    '''
        使用B样本进行插值（B样条插值点：0-0.5-1）
    :param ordinates: 原机翼坐标
    :param fx: 横坐标分布
    :param N: 上/下表面翼型点数
    :return:
    '''

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    np.set_printoptions(threshold=np.inf)  # 打印不省略

    # 定义N为上表面的数据点个数
    N = 65

    # get x, y
    x = ordinates[:, 0].flatten()
    y = ordinates[:, 1].flatten()

    # print("ordinates")
    # print(ordinates)
    # for it in ordinates:
    #     print('''[{},{}],'''.format(it[0], it[1]))

    # B样条插值
    tck, u = spi.splprep([x, y], s=0, per=0)  # 注：per==0不封闭，per==1封闭
    xi, yi = spi.splev(fx, tck)

    ordinates_bspline = np.empty(shape=(len(xi), 2))
    ordinates_bspline[:, 0] = xi
    ordinates_bspline[:, 1] = yi

    # 再次插值得到横坐标相同的分布
    # 插值方法：三次样条插值
    x_up = produce_x_v7(N=N)
    x_low = produce_x_v7(N=N)
    ordinates_CbSp = np.empty(shape=(N * 2 - 1, 2))

    # 处理B样条插值后的数据
    max_id = np.argmax(ordinates_bspline[:, 0])
    ordinates_bspline_tmp = np.concatenate([ordinates_bspline[max_id::],
                                            ordinates_bspline[0:max_id]], axis=0)

    # 剔除异常点
    error_id = []
    for i in range(ordinates_bspline_tmp.shape[0]):
        if np.abs(ordinates_bspline_tmp[i, 0] - 1) < 1e-7:
            error_id.append(i)
    ordinates_bspline_tmp = np.delete(ordinates_bspline_tmp, error_id, axis=0)

    # 将翼型规范化
    ordinates_bspline_tmp = translation(ordinates_bspline_tmp)

    # 翼型数据截断
    ordinates_bspline_tmp = set0(ordinates_bspline_tmp)
    # print("tmp")
    # print(ordinates_bspline_tmp)

    # 剔除重复值，避免三次样条插值出现nan的情况
    equal_list = []
    for i in range(ordinates_bspline_tmp.shape[0] - 1):
        if ordinates_bspline_tmp[i, 0] == ordinates_bspline_tmp[i + 1, 0]:
            equal_list.append(i)
    ordinates_bspline_tmp = np.delete(ordinates_bspline_tmp, equal_list, axis=0)

    # 前缘顶点
    key = np.argmin(ordinates_bspline_tmp[:, 0])
    # print(ordinates_bspline_tmp)
    # 上表面坐标点
    ordinates_CbSp_up = ordinates_bspline_tmp[0:key + 1]
    ordinates_CbSp_up = ordinates_CbSp_up[::-1]

    # 下表面坐标点
    ordinates_CbSp_low = ordinates_bspline_tmp[key::]

    # pd.DataFrame(np.reshape(ordinates_CbSp_up, newshape=(-1, 2))).to_csv("test_bspline_up.csv", index=None, columns=None)
    # pd.DataFrame(np.reshape(ordinates_CbSp_low, newshape=(-1, 2))).to_csv("test_bspline_low.csv", index=None, columns=None)
    # pd.DataFrame(np.reshape(ordinates_bspline, newshape=(-1, 2))).to_csv("test_bspline.csv", index=None, columns=None)

    # test up
    # for i in range(ordinates_CbSp_up.shape[0] - 1):
    #     if ordinates_CbSp_up[i, 0] > ordinates_CbSp_up[i + 1, 0]:
    #         print("up", i, ordinates_CbSp_up[i])
    #
    # # test low
    # for i in range(ordinates_CbSp_low.shape[0] - 1):
    #     if ordinates_CbSp_low[i, 0] > ordinates_CbSp_low[i + 1, 0]:
    #         print("low", i, ordinates_CbSp_low[i])

    pd.DataFrame(np.reshape(ordinates_CbSp_up, newshape=(-1, 2))).to_csv("test/test_up.csv", index=None,
                                                                         columns=None)

    # 上表面
    ipo3_x = spi.splrep(ordinates_CbSp_up[:, 0], ordinates_CbSp_up[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[0:N, 1] = spi.splev(x_up, ipo3_x)  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[0:N, 0] = x_up

    # 下表面
    ipo3_x = spi.splrep(ordinates_CbSp_low[:, 0], ordinates_CbSp_low[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[N - 1::, 1] = spi.splev(x_low, ipo3_x)  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[N - 1::, 0] = x_low

    # 翼型数据截断
    ordinates_CbSp = set0(ordinates_CbSp)

    # show
    # plt.scatter(ordinates[:, 0], ordinates[:, 1], c='b', s=20, label="原机翼")
    # plt.scatter(ordinates_bspline[:, 0], ordinates_bspline[:, 1], c='r', s=5, label="B-样条拟合的机翼")
    # plt.scatter(ordinates_bspline_tmp[:, 0], ordinates_bspline_tmp[:, 1], c='yellow', s=5, label="B-样条拟合的机翼（规范化）")
    # plt.scatter(ordinates_CbSp[:, 0], ordinates_CbSp[:, 1], c='black', s=5, label="三次样条插值拟合的机翼")
    # plt.plot([0, 1], [0, 0], c='orange', linestyle="--")  # 横轴
    # plt.plot([0, 0], [-0.3, 0.3], c='orange', linestyle="--")  # 纵轴
    # plt.ylim(-0.3, 0.3)
    # plt.ylabel('指数')
    # plt.title('机翼数据B样条插值拟合结果')
    # plt.legend()
    # plt.show()

    return ordinates_CbSp


def B_Spline_Interpolation_v3(ordinates, fx, n=65):
    '''
        使用B样本进行插值（B样条插值点：0-0.5-1）
    :param ordinates: 原机翼坐标
    :param fx: 横坐标分布
    :param n: 上/下表面翼型点数
    :return:
    '''

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    np.set_printoptions(threshold=np.inf)  # 打印不省略

    # get x, y
    x = ordinates[:, 0].flatten()
    y = ordinates[:, 1].flatten()

    # B样条插值
    tck, u = spi.splprep([x, y], s=0, per=0)  # 注：per==0不封闭，per==1封闭
    xi, yi = spi.splev(fx, tck)

    ordinates_bspline = np.empty(shape=(len(xi), 2))
    ordinates_bspline[:, 0] = xi
    ordinates_bspline[:, 1] = yi

    ############################################################
    # 将翼型规范化
    ordinates_bspline_tmp = translation(ordinates_bspline)

    # 翼型数据截断
    ordinates_bspline_tmp = set0(ordinates_bspline_tmp)

    # 剔除重复值，避免三次样条插值出现nan的情况
    equal_list = []
    for i in range(ordinates_bspline_tmp.shape[0] - 1):
        if ordinates_bspline_tmp[i, 0] == ordinates_bspline_tmp[i + 1, 0]:
            equal_list.append(i)
    ordinates_bspline_tmp = np.delete(ordinates_bspline_tmp, equal_list, axis=0)

    ############################################################
    # 获取上下表面
    # 前缘顶点
    key = np.argmin(ordinates_bspline_tmp[:, 0])

    # 上表面坐标点
    ordinates_CbSp_up = ordinates_bspline_tmp[0:key + 1]
    ordinates_CbSp_up = ordinates_CbSp_up[::-1]

    # 下表面坐标点
    ordinates_CbSp_low = ordinates_bspline_tmp[key::]
    print("02:", ordinates_CbSp_low[0:2])
    print("tmp", ordinates_CbSp_low.shape)
    print("tmp", ordinates_CbSp_up.shape)

    # print("shape:")
    # print(ordinates_CbSp_low.shape)
    # print(ordinates_CbSp_up.shape)

    # 三次样条插值
    x_up = produce_x_v7(N=n)
    x_low = produce_x_v7(N=n)
    ordinates_CbSp = np.empty(shape=(n * 2 - 1, 2))

    # 判断单调性
    # print("test Monotonicity:", judge_Monotonicity(ordinates_CbSp_low))

    # 上表面
    ipo3_x = spi.splrep(ordinates_CbSp_up[:, 0], ordinates_CbSp_up[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[0:n, 1] = spi.splev(x_up, ipo3_x)[::-1]  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[0:n, 0] = x_up[::-1]

    # 下表面
    ipo3_x = spi.splrep(ordinates_CbSp_low[:, 0], ordinates_CbSp_low[:, 1], k=3)  # 样本点导入，生成参数
    ordinates_CbSp[n - 1::, 1] = spi.splev(x_low, ipo3_x)  # 根据观测点和样条参数，生成插值
    ordinates_CbSp[n - 1::, 0] = x_low

    # show
    # plt.scatter(ordinates[:, 0], ordinates[:, 1], c='black', s=20, label="原机翼")
    # plt.scatter(ordinates_bspline[:, 0], ordinates_bspline[:, 1], c='purple', s=5, label="B-样条拟合的机翼（未规范化）")
    # plt.scatter(ordinates_bspline_tmp[:, 0], ordinates_bspline_tmp[:, 1], c='pink', s=5, label="B-样条拟合的机翼（规范化）")
    # plt.scatter(ordinates_CbSp_up[:, 0], ordinates_CbSp_up[:, 1], c='red', s=5, label="B-样条拟合的机翼（上）")
    # plt.scatter(ordinates_CbSp_low[:, 0], ordinates_CbSp_low[:, 1], c='blue', s=5, label="B-样条拟合的机翼（下）")
    # plt.scatter(ordinates_CbSp[:, 0], ordinates_CbSp[:, 1], c='black', s=5, label="三次样条插值拟合的机翼")
    # plt.plot([0, 1], [0, 0], c='orange', linestyle="--")  # 横轴
    # plt.plot([0, 0], [-0.3, 0.3], c='orange', linestyle="--")  # 纵轴
    # plt.ylim(-0.3, 0.3)
    # plt.ylabel('指数')
    # plt.title('机翼数据B样条插值拟合结果')
    # plt.legend()
    # plt.show()

    return ordinates_CbSp
